{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FiPy 1D two-phase flow in porous mediaq, 11 October, 2019\n",
    "Different approaches:\n",
    "  * Coupled\n",
    "  * Sequential\n",
    "  * ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fipy import *\n",
    "\n",
    "# relperm parameters\n",
    "swc = 0.1\n",
    "sor = 0.1\n",
    "krw0 = 0.3\n",
    "kro0 = 1.0\n",
    "nw = 2.0\n",
    "no = 2.0\n",
    "\n",
    "# domain and boundaries\n",
    "k = 1e-12 # m^2\n",
    "phi = 0.4\n",
    "u = 1.e-5\n",
    "p0 = 100e5 # Pa\n",
    "Lx = 100.\n",
    "Ly = 10.\n",
    "nx = 100\n",
    "ny = 10\n",
    "dx = Lx/nx\n",
    "dy = Ly/ny\n",
    "\n",
    "# fluid properties\n",
    "muo = 0.002\n",
    "muw = 0.001\n",
    "\n",
    "# define the fractional flow functions\n",
    "def krw(sw):\n",
    "    res = krw0*((sw-swc)/(1-swc-sor))**nw\n",
    "    return res\n",
    "\n",
    "def dkrw(sw):\n",
    "    res = krw0*nw/(1-swc-sor)*((sw-swc)/(1-swc-sor))**(nw-1)\n",
    "    return res\n",
    "\n",
    "\n",
    "def kro(sw):\n",
    "    res = kro0*((1-sw-sor)/(1-swc-sor))**no\n",
    "    return res\n",
    "\n",
    "def dkro(sw):\n",
    "    res = -kro0*no/(1-swc-sor)*((1-sw-sor)/(1-swc-sor))**(no-1)\n",
    "    return res\n",
    "\n",
    "def fw(sw):\n",
    "    res = krw(sw)/muw/(krw(sw)/muw+kro(sw)/muo)\n",
    "    return res\n",
    "\n",
    "def dfw(sw):\n",
    "    res = (dkrw(sw)/muw*kro(sw)/muo-krw(sw)/muw*dkro(sw)/muo)/(krw(sw)/muw+kro(sw)/muo)**2\n",
    "    return res\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "sw_plot = np.linspace(swc, 1-sor, 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the relative permeability and fractional flow curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3yV5f3/8dcne+9Fwggj7CEQUBFURAVxIIoVd9VKtbXDb4d+tbZ+21+tbb/229o6iqNqHRT3AATrRClK2ENG2CEJZJO9zvX74zrBCIGcQJL7nJPP8/E4j7Nuct6G+ObOdd/3dYkxBqWUUr4vwOkASimlOocWulJK+QktdKWU8hNa6Eop5Se00JVSyk8EOfXBSUlJJjMz06mPV0opn7R69epiY0xyW+85VuiZmZnk5OQ49fFKKeWTRGTv8d7TIRellPITWuhKKeUntNCVUspPaKErpZSf0EJXSik/0W6hi8gzInJIRDYd530RkUdEJFdENojIuM6PqZRSqj2e7KE/C8w4wfsXAVnu2zzg8VOPpZRSqqPaLXRjzKdA6Qk2mQU8b6yVQJyI9OqsgMeoPAhL7oGmhi77CKWU8kWdMYaeAexv9TzP/doxRGSeiOSISE5RUdHJfdr+lfDF4/D+/Sf355VSyk91RqFLG6+1uWqGMWa+MSbbGJOdnNzmlavtGz4LzvgefPEEbHz15L6GUkr5oc4o9DygT6vnvYH8Tvi6x3fBr6HvmfD2D+Dgli79KKWU8hWdUehvAze6z3Y5A6gwxhR0wtc9vsBguOpZCImChTdA3eEu/TillPIFnpy2+DLwH2CIiOSJyK0icruI3O7eZDGwC8gFngS+12VpW4tOs6Veuhve+h7o2qhKqR6u3dkWjTHXtPO+Ab7faYk6IvMsO/yy7D5Y8Qic9SNHYiillDfw/StFz/y+PVD67wdg93Kn0yillGN8v9BFYNajkDgIXr0ZDnft8VillPJWvl/oAKHRcPUL0FADC2+EpnqnEymlVLfzj0IHSB4Csx+HvFWw5G6n0yilVLfzn0IHO5Y++S5Y/Q9Y/ZzTaZRSqlv5V6EDnHc/DDwPFv8U8nTNUqVUz+F/hR4QCFc+DdG94F83QNUhpxMppVS38L9CB4hIgLkvQm0ZLLwJmhudTqSUUl3OPwsdIG0UXPZX2LcClv3C6TRKKdXl2r1S1KeNvgry18LKRyF9LIyZ63QipZTqMv67h97igl9D5hR4+4dwYI3TaZRSqsv4f6EHBsFVz0FUKiy4zq54pJRSfsj/Cx0gMhGueQnqyu10u3olqVLKD/WMQgd7kPTyx2D/F/YcdZ1uVynlZ/z7oOjRRsyGwk2w/H8hbTRMvM3pREop1Wl6zh56i6n3weAZ8N49Ot2uUsqv9LxCDwiAK+ZDwgB45SYo2+t0IqWU6hQ9r9ABwmJh7svQ3AQLroX6KqcTKaXUKeuZhQ6QNAjmPAOHtsAb3wWXy+lESil1SnpuoQNknQ8X/ha2vgsfP+h0GqWUOiU96yyXtpxxh91L//SPkDwURs1xOpFSSp2Unr2HDnZN0ov/BH0nwVvfh7zVTidSSqmTooUOEBQCV/8TolLsQVJdaFop5YO00FtEJsE1C6ChypZ6Q43TiZRSqkO00FtLHQFXPAn56+zwi04PoJTyIVroRxs6E87/FWx+HT7+ndNplFLKY3qWS1vO+jEU58Inv4fEQTD6W04nUkqpdukeeltE4JL/g36T7dDLvpVOJ1JKqXZpoR9Py5kvsX3sQdLS3U4nUkqpE9JCP5GIBLh2Ibia4aWrobbc6URKKXVcWujtSRpk99RLd8Ir34bmRqcTKaVUmzwqdBGZISLbRCRXRO5p4/1YEXlHRNaLyGYRubnzozqo/9lwyZ9h10ew5Od6OqNSyiu1e5aLiAQCjwIXAHnAKhF52xizpdVm3we2GGMuFZFkYJuIvGiMaeiS1E4YdwOU7IDP/2LnUp/0A6cTKaXUN3iyhz4RyDXG7HIX9AJg1lHbGCBaRASIAkqBpk5N6g2mPQDDZ8Gy+2HLW06nUUqpb/Ck0DOA/a2e57lfa+1vwDAgH9gI/MgYc8wE4yIyT0RyRCSnqKjoJCM7KCAAZv8demfD6/Ng/yqnEyml1BGeFLq08drRg8jTgXVAOnAa8DcRiTnmDxkz3xiTbYzJTk5O7nBYrxAcbud8iU6Dl+fq6YxKKa/hSaHnAX1aPe+N3RNv7WbgdWPlAruBoZ0T0QtFJsF1r4FphhevgppSpxMppZRHhb4KyBKR/iISAswF3j5qm33ANAARSQWGALs6M6jXSRoEc1+C8r3wr+uhqd7pREqpHq7dQjfGNAF3AkuBr4CFxpjNInK7iNzu3uw3wCQR2Qh8ANxtjCnuqtBeo98kuPxx2Pu5zs6olHKcR5NzGWMWA4uPeu2JVo/zgQs7N5qPGDUHyvfBB/8DcX1h2i+dTqSU6qF0tsXOMPkuO/Sy/GGISYcJ33E6kVKqB9JC7wwiMPNhqDwIi38G0b1g6MVOp1JK9TA6l0tnCQyCOU9D+lh49RbY/6XTiZRSPYwWemcKibSzM8ak29kZi3OdTqSU6kG00DtbZBJc/xpIALxwhR2GUUqpbqCF3hUSBsB1C6G6CF66CuornU6klOoBtNC7SsZ4uOo5KNwEC2+EJv+ZeFIp5Z200LvS4Avhskdg54fw5h3gOma+MqWU6jR62mJXG3s9VB2yFx5FJsOM39nTHJVSqpNpoXeHyXfZUv/icYhKgSn/5XQipZQf0kLvDiIw/UGoKf56T33cDU6nUkr5GS307hIQALMes1PtvvNDiEiEoTOdTqWU8iN6ULQ7BYXAt56HXqfBqzfD3hVOJ1JK+REt9O4WGgXXvQKxveGluVC40elESik/oYXuhMgkuOFNW+7/vAJKdjqdSCnlB7TQnRLXx5a6aYbnL4eKA04nUkr5OC10JyUPtvO+1JbBP2dDdYnTiZRSPkwL3WnpY+HaBXaBjBev1HlflFInTQvdG2ROtvO+FGyAl6+BxjqnEymlfJAWurcYMgNmPwF7PrMLZDQ3Op1IKeVjtNC9yehvwcw/wrZF7sm8mp1OpJTyIXqlqLeZeBs0VMG/H4DgcLj0EZ3MSynlES10bzT5LmiogU//AMERMOMhLXWlVLu00L3V1HuhoRpWPmrXKp32S6cTKaW8nBa6txKB6b+FxhpY/rDdUz/7p06nUkp5MS10byYCF//JlvqHv7F76mfc4XQqpZSX0kL3di3T7jbWwHv3QFAoZN/idCqllBfS0xZ9QWAQXPkMZE2Hd++CNf90OpFSygtpofuKlrnUB54Hb/8A1r3sdCKllJfRQvclwWEw9yXofza89T3Y8IrTiZRSXsSjQheRGSKyTURyReSe42xzroisE5HNIvJJ58ZURwSHwzULoO8keOO7sPkNpxMppbxEu4UuIoHAo8BFwHDgGhEZftQ2ccBjwGXGmBHAVV2QVbUIiYBr/wW9J8Br34Gv3nU6kVLKC3iyhz4RyDXG7DLGNAALgFlHbXMt8LoxZh+AMeZQ58ZUx2hZyi59LLzybdi2xOlESimHeVLoGcD+Vs/z3K+1NhiIF5GPRWS1iNzY1hcSkXkikiMiOUVFRSeXWH0tLMYukJE2Cv51A2xd7HQipZSDPCn0tiYRMUc9DwLGAxcD04H7RWTwMX/ImPnGmGxjTHZycnKHw6o2hMXCDW/YUl94I2xd5HQipZRDPCn0PKBPq+e9gfw2tnnPGFNtjCkGPgXGdE5E1a7wOLjxTeg12pb6V+84nUgp5QBPCn0VkCUi/UUkBJgLvH3UNm8BU0QkSEQigNOBrzo3qjqhlj31XqfZMfUtR/8VKaX8XbuFboxpAu4ElmJLeqExZrOI3C4it7u3+Qp4D9gAfAk8ZYzZ1HWxVZtaSj19HLx6M2x5y+lESqluJMYcPRzePbKzs01OTo4jn+336g7Di3MgLwfmPA0jZjudSCnVSURktTEmu6339EpRf9Ry9kvvCfDqrXpFqVI9hBa6vwqNtqXebxK8fptO6KVUD6CF7s9Co+DahTBwKrx9J6x6yulESqkupIXu70IiYO7LMPgiWPQT+M+jTidSSnURLfSeIDjMTr077DJYeq9d0k4p5Xe00HuKoBCY8w8YdRV88Gv48Lfg0BlOSqmuoUvQ9SSBQTD773YZu0//AA1VMP1Bu3apUsrnaaH3NAGBcOlfISQKVj4G9Yfh0kfs60opn6aF3hMFBMCMhyA0xu6p11fBFU/aYRmllM/SQu+pROC8++xFSMt+YYdfvvVPe1aMUson6UHRnm7SD+yQS+4H8MKVUFfhdCKl1EnSQlcw/iY750vel/DcpVCli48o5Yu00JU18kqY+xIUbYN/zICyvU4nUkp1kBa6+trg6XDDm1BdBM9Mh4NbnE6klOoALXT1Tf3OhJuX2IuO/jED9n3hdCKllIe00NWxUkfArcsgIgmenwXblzmdSCnlAS101bb4fnDLUkgeDC/PhfULnE6klGqHFro6vqhkuOldyDwL3vgufP6Izv+ilBfTQlcnFhYD171ql7F7/3547x5wNTudSinVBr1SVLUvKBSufAai02Hlo1BZALPn22l5lVJeQwtdeSYgAGY8CDHpsOw+e/HR3BchIsHpZEopNx1yUR0z6U6Y8wwcyIFnZkD5fqcTKaXctNBVx428Eq5/HSoL4anzoXCj04mUUmihq5PVfwrc8p6dR/2ZGbDjfacTKdXjaaGrk5c6HL7zb0gYAC9dDauedjqRUj2aFro6NTHpdqqAQefDov+CpfeBy+V0KqV6JC10depCo+xMjRNug//8DRbeAA01TqdSqsfRQledIzAIZv4Rpv8Oti6C5y6BqkNOp1KqR9FCV51HBM78Hlz9gp1698nzoHCT06mU6jG00FXnG3YJ3LwYXE12XvVtS5xOpFSPoIWuukbGOLjtQ0gcBC9fA5//RSf2UqqLeVToIjJDRLaJSK6I3HOC7SaISLOIzOm8iMpntZwBM3wWvP9LeOv70FTvdCql/Fa7hS4igcCjwEXAcOAaERl+nO1+Dyzt7JDKh4VEwJx/wDn3wLoX7YIZ1cVOp1LKL3myhz4RyDXG7DLGNAALgFltbPcD4DVAT21Q3xQQAFP/G658GvLXwpNTdboApbqAJ4WeAbSegSnP/doRIpIBzAaeONEXEpF5IpIjIjlFRUUdzap83ag59mBpcyM8fSFsfsPpREr5FU8KXdp47eijW38G7jbGnHDlA2PMfGNMtjEmOzk52dOMyp9kjId5H0PqSHjl2/DBr3XBDKU6iSfzoecBfVo97w3kH7VNNrBARACSgJki0mSMebNTUir/Ep0G334XFv8Mlj9sz1W/8kkIi3U6mVI+zZM99FVAloj0F5EQYC7wdusNjDH9jTGZxphM4FXge1rm6oSCQuHSv8DFD8POD+DJaVC8w+lUSvm0dgvdGNME3Ik9e+UrYKExZrOI3C4it3d1QOXHRGDCd+DGt6G2zF5ZunWx06mU8lliHLrYIzs72+Tk5Djy2coLle+3k3rlr4UpP4Wp99q51pVS3yAiq40x2W29p1eKKu8Q1wdufg/G3QjL/xdeuBKqS5xOpZRP0UJX3iM4DC77K1z6COxdAfPPgQNrnE6llM/QQlfeZ/xNdnk7sJN7rX7O2TxK+QgtdOWdMsbBvE+g31nwzg/hze/pohlKtUMLXXmvyES4/jU4++ew7iV4ahoUbXc6lVJeSwtdebeAQDjvPlvsVYdg/rmwYaHTqZTySlroyjcMmga3L4deY+D12+CdH0FjrdOplPIqWujKd8Skw03vwOS7YPWz8NQFULLT6VRKeQ0tdOVbAoPg/Afg2lfgcB78/WxYv8DpVEp5BS105ZsGXwi3f2aHYN74Lrw+D+ornU6llKO00JXviu1th2DO/W/Y+IrdW89f63QqpRyjha58W0AgnHsPfHuRXa/0qQtgxV/B5XI6mVLdTgtd+Yd+k+wQzODpsOwX8OIcqCx0OpVSxzh4uI5DlXVd8rW10JX/iEiAq1+wc6zv/RweOxO+esfpVErR2Oxi6eZCbn12FZMe+pCnl+/uks/xZMUipXxHyxzrmVPs+er/uh7G3gAzHoLQKKfTqR4m91AVC3P28/qaPIqrGkiJDuW7Zw/gW9l92v/DJ0ELXfmn5CFw67/h49/BZ/8Hez6DK+ZDn4lOJ1N+rrKukUUbCnhldR6r95YRFCCcNzSFqyf04ZzByQQFdt3AiBa68l9BIXD+ryDrAnj9u3bmxrN/Zm+BwU6nU37E5TKs3FXCK6vzWLKpgLpGFwOTI7l35lBmj+1NcnRot+TQQlf+r98kuOMzWHI3fPJ72P4eXP4EpA53OpnycftKanh1TR6vrc7jQHkt0WFBXDGuN1eN781pfeIQkW7No4WueoawWJj9BAy9GN75sV0849x7YNKP7NWnSnnocF0jizcU8NqaPFbtKUMEJg9K4uczhjB9RBphwc4tnag/yapnGXYp9D0TFv0XfPBruyj15Y9D8mCnkykv1uwyLN9RxGtrDrBscyH1TS4GJEfys+lDmD02g/S4cKcjAlroqieKTIKrnoNNr8Hin8Lfp8B598MZd+jC1OoIYwyb8w/z5toDvLU+n6LKemLDg/lWdh+uGJfhyJBKe7TQVc8kAqPmQOZkOwSz7D7Y8hbM+ps9Q0b1WPnltby1Lp831uax/WAVwYHCuUNSuGJsBucNSyE0yHv/0RdjjCMfnJ2dbXJychz5bKW+wRi7aMZ7d0NDNZzzczjrx3omTA9SUdPIkk0FvLUun5W7SzAGxvWNY/a43lwyqhfxkSFORzxCRFYbY7Lbek/30JUSgTFXw8DzYMnP4MP/B5vftHvr6WOdTqe6SF1jMx9uPcSbaw/w8bYiGppd9E+K5EfTspg9NoN+iZFOR+wwLXSlWkQlw1XPwsg5sOgn8OQ0mHSnnc0x2DsOeqlT09TsYsXOEt5en8/STYVU1jeRHB3K9Wf04/Kx6YzKiPW6cfGO0EJX6mjDLrFj6+/fD5//xY6tX/wnuwye8jkulyFnbxnvrM9n8cYCSqobiAoNYsbINC4/LYMzByYSGOC7Jd6aFrpSbQmPg8v+CqOugnfvgheusHvu0x+E6FSn06l2GGPYeKCCdzcU8O76fPIr6ggLDmDasFQuHZ3OuUOSHT1fvKtooSt1Iv3PhjtW2Plglj8MO9630wmMvxkCdLJSb9JymuGijQUs2lDAvtIaggOFs7OSufuioUwblkpUqH9Xnp7lopSninfYvfU9y6H3BLjkz5A20ulUPZoxhq2FlSzaUMCijQXsLq4mMEA4a1ASl4zuxfThacRG+NfZSnqWi1KdISnLLnm34V+w9F675N3E2+xB0/A4p9P1GC174os3FrBkUyG7i6sJEJg0MIl5Zw9g+og0ErzoNMPu5FGhi8gM4C9AIPCUMeaho96/Drjb/bQKuMMYs74zgyrlFURgzFzIutCe3vjF3+0Vp+f/D4y5RodhukjLmPjijYUs2VTA3pIaAgOEMwckctuUAVw4IpWkqO6Z0dCbtTvkIiKBwHbgAiAPWAVcY4zZ0mqbScBXxpgyEbkIeMAYc/qJvq4OuSi/kL/OTh+Qtwp6T4SZf4T005xO5ReaXYbVe8tYsqmApZsKya+oIyhAmDQoiZkj07iwh+6Jn+qQy0Qg1xizy/3FFgCzgCOFboxZ0Wr7lUDvk4+rlA9JPw1uWQbrX4L3fwXzz4Xsm2HqLyAy0el0PqehycXKXSW8t7mQZZsPUlxVT0hQAGdnJXHXBYO5YHgqcRE9r8Q95UmhZwD7Wz3PA060930rsKStN0RkHjAPoG/fvh5GVMrLBQTA2Oth6CXw0YOw6inY+BqcezdMuM0utKGOq7q+iU+2F7F0cyEfbj1EZV0TESGBTB2awowRaUwdmuL3Z6d0Fk++S22dcd/mOI2ITMUW+uS23jfGzAfmgx1y8TCjUr4hPA5m/gGyb7EHTZfeC6uehum/hcEz7Pi7AqCosp4Ptx5k6eaDfJZbTEOTi/iIYGaMSGP6iDQmZyX55XniXc2TQs8DWq9o2hvIP3ojERkNPAVcZIwp6Zx4SvmglKFww+v2nPWl98LLc6H/OfaipB56mqMxhp1FVSzbcpB/bznI2v3lGAMZceFcf3o/LhyRSna/+C5db7Mn8KTQVwFZItIfOADMBa5tvYGI9AVeB24wxmzv9JRK+aKsC2DAuZDzD/j4QTvv+phrYeq9EJvhdLou19jsYvXeMj746iDvbznInpIaAEZlxHLX+YOZNiyF4b1ifHruFG/TbqEbY5pE5E5gKfa0xWeMMZtF5Hb3+08AvwQSgcfcfzlNxzsKq1SPEhgMp8+zc68vfxi+nA+bXoXTb4fJd/nd+etl1Q18sr2ID7Ye4pNthzhc10RIYABnDkzkO1MGMG1YCr1idaKzrqJXiirVncr22vPXNy6E8HiY8lN7cVKQb55DbYxh28FKPtpaxEdbD5GztxSXgaSoEKYOSWHasFQmZyXpQc1OdKLTFrXQlXJCwXp7muOujyC2r12wevTVPrFgdU1DE5/nlvDRtkN8vPUQ+RV1AAzrFcP5w2yJj86IJcBPZjD0NlroSnmrnR/Cvx+wBZ+YBVP/G4bP9qorTlsOaH68rYhPthfxxa5SGppdRIYEMjkrialDUjh3SAppsWFOR+0RtNCV8mbGwFfvwEe/haKtkDoSpt4HQy5y7FTHyrpGVuws4ZPtRXyyrYgD5bUADEqJ4pzByZw3NIUJmQmEBHnPPzw9hRa6Ur7A1WznhfnoQSjbDRnj7RkxA6d1ebE3u+xcKcu3F7F8RzFr9pXR5DJEhgRy1qAkzh2SwtmDk+gdH9GlOVT7tNCV8iXNjbD+ZfjkD1CxH9LHwTl3w+DpnVrsB8pr+WxHEZ/uKObz3GLKaxoBGJkRw9lZyUzJSmZ8v3jdC/cyOn2uUr4kMBjG3Qij59o5Ypb/CV6+GtJGwzk/hyEXn9QYe0VNI//ZVcxnucV8nlvC7uJqAFJjQjl/WCpTspKYPCiJRJ210GfpHrpS3q65ETYstOexl+6ElOEw5Scw/PITnhVT29BMzt5SVuwsYUVuMRsPVOAyEBESyBkDEjlrkC3wwalRenGPD9EhF6X8QXMTbH4DPv0jFG+D+Ew480447ToIiaChycX6vHJW5JawYmcxa/eV09DsIihAGNMnjsmDkpiclcSY3nE6jOLDtNCV8icuF2xbjOuzPxNwYBW1wXEsCr+Mh8vOpqAxAhEYkR7DpIFJnDkwkQmZCXphjx/RMXSl/EB9UzMb8ir4cncpK3clk7PvJ4xq2sx3m99lTuPzXBa8kPwhV5Ew9QfEZAxxOq5ygBa6Ul6qpqGJNXvL+XJ3CV/sLmXt/nIamlwADEmN5uoJfTljwGmM7f9DqNpByOePkLnpZch9AYbMhDPugMzJOm1vD6JDLkp5iaLKelbvLWXVnjJy9pSyKf8wzS5DgMDIjFgmZiYwsX8CEzITiD/e0muHC+wCGznPQG0ppI6yxT7ySgjWKzn9gY6hK+VlXC7DruIqVu8tI2dPGTl7y46cRhgaFMCYPnFk94tnYv8ExveLJzosuGMf0Fhrz4xZ+TgUfQWRyTD2Brs8XpyuFubLtNCVclh1fRPr88pZs7eM1XvLWLOvnIpaeyFPfEQw4/slMCEznuzMBEZmxBAa1Emr9RgDuz6GL/4OO5ba54OnQ/atMGgaBOiqQL5GD4oq1Y1a9r7X7Ctn7b5y1u4rY/vBSlzufafBqVHMHJXGuL7xjO8XT/+kyK47D1wEBk61t/L9sPpZWPM8bH/P7qmPv9muhxqV0jWfr7qV7qErdYoOVdaxfn8F6/eXsz6vnHX7y6msawIgOiyI0/rEMbZvPGP7xjGuTzyxER0cPulsTQ2w9V07zr5nOQQEQdZ0GHcDDLrAJ6bw7cl0D12pTlJR08im/Ao25FWwIa+c9fvLj8wHHhggDEmN5pLR6Yzra0t8QFKk980LHhQCI6+wt6LtsPafdu6YbYsgKhXGXGP32pOynE6qOkj30JU6jsN1jWw+cJhNByrYcKCCjXnlR9bFBOibEMFpfeIY0yeOMb1jGZEeS3iIj45JNzfCjmWw9gXYvhRMM/SeCKO/BSNmQ2SS0wmVmx4UVaodpdUNbDpQwab8Clvi+RXsbVXeGXHhjMqIZVTvWEb3jmVURixxEcc5ddDXVRbC+gX2LJlDm+2QzMBpttyHzIQQnULXSVroSrm5XIZ9pTVsKTjMlvzDR+4LD9cd2aZPQjgj02MZmRHL8PQYRmXEktRTZyAs3GTXP934Khw+ACFRMPRiOzHYwPP03HYHaKGrHqmitpFthZVsLTzM1sJKthYcZlthJdUNzYAd8x6UHMXw9BiG9YpmZEYsI3rFOn/Q0hu5XLD3c9jwL7u6Ul05hETbVZVGXG734LXcu4UWuvJrdY3N5B6qYvvBSrYftPfbCiuPLJsGEBsezJC0aIalRTM8PYbhvWLJSo0iLNhHx7yd1NwIuz6BLW/as2Vqy+ye++AZMHQmDDofwmKdTum3tNCVX6htaGZnURU7i6rYcbCKHYcq2XGwij0l1UfO8Q4OFAYmRzEkLZqhaTEMTYtmaK9o0mLCdM7vrtDcCLs/dZf7IqgpsWPumZPtePvgGRDfz+mUfkULXfkMYwyl1Q3sLKpml7u8dxZVk3uoiv1lNbT8uAYGCP0SIxicEs3gtGiGpEYzJC2KfomRBAfqXN+OcDVDXg5sWwzbltg52wFSRkDW+XZYpu8ZENRDj0d0Ei105XVqGprYU1zDnpJqdhdXs6uomt3FVewqrj6ytiXYeU36J0UyKCWKrJRoe58aRWZipC7S4O1Kdtpi3/4e7FsJrkYIjrR774Om2YJPHKizQXaQFrpyRFV9E3tLqtlbUuO+VbPH/bygou4b26bGhJKZGMnAlCgGJkcxMDmSgclRpMeFE+htF+aojquvhN3LYecHkPsBlO22r8f0tgXffwpkTtHhGQ/olaKqSzQ2uyisqGN/aQ37SmvYX1bDvtJa9pfWsL+0hpLqhm9snxQVQr/ESM4ckEj/pEj6J0fSPymSzMRIInVFHf8WGh0LiMYAAAnVSURBVG0PmA6daZ+X7ISdH9qpB3Lfhw0L7OuxfW259z0T+pwOiYNOakHsnkr30NVxVdc3UVBRS355HfnltRworyWvrJYDZbXkldVQeLjuyMFIsOPaGXHh9EkIp29CBH0SIshMjKRfYgT9EiN1GTTVNmOgaKvdg9+zHPZ8ZudyBwiPh94ToM9Ee+VqxngIjXI2r8N0D10do6q+icKKWgor6imoqOXg4ToKKuoorKgjv8IWeMv0ri0CBHrFhpMRH84ZAxLpHW8f94m35d0rNowgPSCpOkoEUobZ2+nz7DnvJbmw/wvI+xL2f2mnJbAb2zlmep0G6afZ+16j7W8ASvfQ/YnLZSiraaC4qoHiqnqKKu3t4OE6DrnvW563XFzTWnxEMKkxYWTEhdMrLoz0uHDSY8PpFWsfp8WG6Rkkyhm1ZfYMmgOrIX8dFKyDygL3mwIJAyB1OKS0uiUM8MuZI3UP3Uc1uwwVtY2UVjdQXtNAabW9lVQ3UFLVQGl1/ZHHxVX2cbPr2H+gw4IDSI0JIyU6lGHpMZwzJJnUmDB6xYZ9414vslFeKzwesi6wtxaVB6FgvS33gvVw6Ct7Lryx664SGAJJg+04fOJASBj49eOIRL88u8ajQheRGcBfgEDgKWPMQ0e9L+73ZwI1wLeNMWs6OatPqm9qprKuyX1r5HCtva+sa+JwXSPlNY1U1H59K69t5LC7xA/XNXK8X6AiQgJJjAohITKUtNgwRmbEkBwdSlJU6Dfuk6NDiQ4N0otqlP+JToXoC2HwhV+/1lgLxdttuR/aYu8LN9rpCkyr30pDYyG+L8T2cd96Q5z7cUy6XbIv0PemgGi30EUkEHgUuADIA1aJyNvGmC2tNrsIyHLfTgced997JZfL0Ohy0dRsaGx20dDkot59a2hy0eB+ra6xmdrGZuoam6lvdFHX1ExtQzM1Dfb1moYmaurt85rGZqrrm6iut+Vd3WAfNzafeEgrQCAmPJjY8GDiwoOJCQ+mT3w4CZEhxEWEEB8R/I3HiVGhJEaG6N60Um0JDodeY+ytteZGKN9nz64pyYXSnXYFp7I99mBsQ+WxXys8wc4PH5Vi7yOTICwOwuPc9/Hux7EQHOG+hUNQmGNn5niyhz4RyDXG7AIQkQXALKB1oc8Cnjd2QH6liMSJSC9jTMGxX+7UfLztEL95dwsGwIDBXl1o7+0whcuYb9y33BpdhqZmF22MSnRYWHAAESFBhAcHEhFib1FhQSRGRhAVGkRUWBCRoUFEhQYRHea+hQa7H9v7mPBgokODvG8BBKX8TWCwHWpJHAhceOz7teVQkQcV++3YfNUh9+2gvd//hZ3WoKHKs88LCrflHhBk120NCAIJ+PrxuJtg0p2d+p8InhV6BrC/1fM8jt37bmubDOAbhS4i84B5AH37ntzK49FhwQxNiwEBsV/TfW+fBwQIgSIEBsg3HgcGCEGBQnBAgL0PDCAowN6HBAUQGnT0fSBhwfY+PCSQsOBAwoICCAsOJDw4UEtYKX8S7t7zTht54u2aG6Guwv4DUFf+9X1THTTUQGONHfZpuXc12aEeV8vN/byL1nD1pNDbaq6j93E92QZjzHxgPtizXDz47GOM72cX1lVKqW4XGGyHXrx0BSdPBnrygD6tnvcG8k9iG6WUUl3Ik0JfBWSJSH8RCQHmAm8ftc3bwI1inQFUdMX4uVJKqeNrd8jFGNMkIncCS7GnLT5jjNksIre7338CWIw9ZTEXe9rizV0XWSmlVFs8Og/dGLMYW9qtX3ui1WMDfL9zoymllOoIvY5bKaX8hBa6Ukr5CS10pZTyE1roSinlJxybPldEioC9J/nHk4DiTozTWbw1F3hvNs3VMZqrY/wxVz9jTHJbbzhW6KdCRHKONx+wk7w1F3hvNs3VMZqrY3paLh1yUUopP6GFrpRSfsJXC32+0wGOw1tzgfdm01wdo7k6pkfl8skxdKWUUsfy1T10pZRSR9FCV0opP+HVhS4iM0Rkm4jkisg9bbw/VET+IyL1IvJTL8p1nYhscN9WiMiYtr6OA7lmuTOtE5EcEZnsDblabTdBRJpFZI435BKRc0Wkwv39Wiciv/SGXK2yrRORzSLyiTfkEpGftfpebXL/XSZ4Qa5YEXlHRNa7v1/dMhusB7niReQN9/+TX4pIO8slecAY45U37FS9O4EBQAiwHhh+1DYpwATgt8BPvSjXJCDe/fgi4AsvyRXF18dNRgNbvSFXq+0+xM7qOccbcgHnAu92x89VB3PFYdf07et+nuINuY7a/lLgQ2/IBdwL/N79OBkoBUK8INcfgV+5Hw8FPjjVz/XmPfQji1MbYxqAlsWpjzDGHDLGrAIavSzXCmNMmfvpSuwKTt6Qq8q4f3qASNpYJtCJXG4/AF4DDnVDpo7k6m6e5LoWeN0Ysw/s/wdekqu1a4CXvSSXAaJFRLA7NaVAkxfkGg58AGCM2QpkikjqqXyoNxf68RaedlpHc90KLOnSRJZHuURktohsBRYBt3hDLhHJAGYDT9B9PP17PNP9q/oSERnhJbkGA/Ei8rGIrBaRG70kFwAiEgHMwP4D7Q25/gYMwy6LuRH4kTHG5QW51gNXAIjIRKAfp7jz582F7tHC0w7wOJeITMUW+t1dmsj9cW281tZC3W8YY4YClwO/6fJUnuX6M3C3Maa5G/K08CTXGuy8GWOAvwJvdnkqz3IFAeOBi4HpwP0iMtgLcrW4FPjcGFPahXlaeJJrOrAOSAdOA/4mIjFekOsh7D/M67C/oa7lFH9z8GjFIod468LTHuUSkdHAU8BFxpgSb8nVwhjzqYgMFJEkY0xXTl7kSa5sYIH9jZgkYKaINBljurJA281ljDnc6vFiEXnMS75feUCxMaYaqBaRT4ExwHaHc7WYS/cMt4BnuW4GHnIPN+aKyG7smPWXTuZy/3zdDOAeDtrtvp28rj5ocQoHFYKAXUB/vj6oMOI42z5A9x0UbTcX0Be7vuokb/p+AYP4+qDoOOBAy3Nv+Ht0b/8s3XNQ1JPvV1qr79dEYJ83fL+wwwcfuLeNADYBI53O5d4uFjtGHdnVf4cd+H49Djzgfpzq/rlP8oJccbgPzgK3Ac+f6ud67R668WBxahFJA3KAGMAlIj/GHkk+fNwv3A25gF8CicBj7r3OJtPFM755mOtK4EYRaQRqgauN+6fJ4VzdzsNcc4A7RKQJ+/2a6w3fL2PMVyLyHrABcAFPGWM2OZ3LvelsYJmxvz10OQ9z/QZ4VkQ2YodC7jZd+1uWp7mGAc+LSDP2rKVbT/Vz9dJ/pZTyE958UFQppVQHaKErpZSf0EJXSik/oYWulFJ+QgtdKaX8hBa6Ukr5CS10pZTyE/8fp83zLvH1rpgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXiU9b3+8fcnIQlbgMSEnQSQRXYIIaCeWq1aUY/FrRVQQVyoVjkuP1s9ntau52q17altXRAtKm6gFautVrRUxVpFEgirLCEIJCBJIBvZJpn5/v5ItCkGmcAkz8zkfl1XruSZeZK5HWZuv/N9NnPOISIikS/G6wAiIhIaKnQRkSihQhcRiRIqdBGRKKFCFxGJEp28euCUlBQ3ePBgrx5eRCQi5eTklDjnUlu6z7NCHzx4MNnZ2V49vIhIRDKz3Ue7T1MuIiJRQoUuIhIlVOgiIlFChS4iEiVU6CIiUeKYhW5mi82syMw2HeV+M7PfmVmemW0ws4zQxxQRkWMJZoT+JDD9S+4/Hxje9DUfeOTEY4mISGsdcz9059wqMxv8JavMAJa4xvPwfmhmvcysn3Nuf4gyioiEtUDAUVFbT2l1PaXVPsqr6ymr8eFrCODzO+obAtT7G798fkdmehJnjGjx2KATEooDiwYAe5stFzTd9oVCN7P5NI7iSUtLC8FDi4i0vboGP7sPVrP3UNNXaQ0FpdXsPVTD/vIaymrqac2lJW786slhW+jWwm0t/qc55xYBiwAyMzN1ZQ0RCTvlNfVs2VfBlv0VbN5XzpZ9FeQVHaYh8K/K6hwXw6CkrgxK7kpGei+Su8bTq2s8vbrGkdT0vWeXODrHxRIXG0N8bAydYo242BjiYg2zlmrzxIWi0AuAQc2WBwL7QvB3RUTaXMnhOj7YeZAP8g/ywc6D7Cqp+vy+3okJjO7fg7NH9WZEn0TSkrsyMKkrKd3j26yUT0QoCv1V4BYzWwpMBco1fy4i4aq23s97O0p4P6+ED3YeZNuBSgC6J3Ri6pBkLp88kDH9ezCmf09SExM8Tts6xyx0M3seOBNIMbMC4IdAHIBzbiHwOnABkAdUA/PaKqyIyPGo8fl5d3sRr238lL9/fIAqn5/OcTFMGZzMjEn9Oe3kFMb270Gn2Mg+NCeYvVxmHeN+B9wcskQiIiFQ1+Bn5cdFvLZxP29vLaLa5ye5WzzfmNif88f2Y+rQZBI6xXodM6Q8O32uiEhb2H2wiuc+2sOL2QUcqvKR0j2eSyYN4MJx/cgakhzxo/Avo0IXkYjX4A+wcmsRz3y4m/d2lBAbY5wzqjdXTk3n9GEpxMaE3wbMtqBCF5GIVVXXwDMf7uaJ9z/h04pa+vbozO3njOCKKYPo27Oz1/HanQpdRCJOZW09Sz7YzePv5VNaXc9pJ5/Ej2eM4exTekf1lMqxqNBFJGKU19Tz5Puf8Id/5FNR28BZI1NZcPZwMtKSvI4WFlToIhL2anx+Hnsvn8dW5VNZ18C5o/uw4GvDGD+wl9fRwooKXUTClnOO1zbu5+evb6WwrIbzxvTh1rNHMLp/D6+jhSUVuoiEpU2F5fzkz1v46JNDjO7Xg//71gSmDj3J61hhTYUuImGl5HAdv35zG0vX7CWpazw/v3Qc38oc1GF2PTwRKnQRCQvOOV5eV8iPXt1Mtc/PdacPYcHZw+nZJc7raBFDhS4iniuqqOWelzfyt4+LmJyexH2XjWdY7+5ex4o4KnQR8UzzUXldQ4DvXziKeacP0fTKcVKhi4gnDlTUcs/yjazc2jgq/+Xl4xmaqlH5iVChi0i7e2vLAe58cT219X6NykNIhS4i7abeH+BXK7bx6Kp8xg7owe9mTtKoPIRU6CLSLvaX13DLc+vI2V3KVdPS+P6Fo+kcF13nI/eaCl1E2ty724u5fVkutfV+fjtzIjMmDvA6UlRSoYtIm/EHHA/8bTsPvp3HiN6JPHRlhnZHbEMqdBFpE4frGrj1+XWs3FrENycP5CczxtIlXlMsbUmFLiIht6+shmufXMOOosP8ZMYY5pw62OtIHYIKXURCav3eMq5fkk2tz8/ia6bw1RGpXkfqMFToIhIyr2/cz+3LcklNTODZ66cyok+i15E6FBW6iJww5xwPv7OTX67YxuT0JB69ejIp3RO8jtXhqNBF5IT4A47v/2kjz3+0lxkT+3PfZeO1f7lHVOgictx8DQFuX5bLaxv3c/NZJ3Pn10dipkP4vaJCF5HjUu1r4MZn1rJqezH/c8EobjhjqNeROjwVuoi0WnlNPdc9uYa1e0q577JxXDElzetIggpdRFqpuLKOuYs/YkdRJQ/OzuCCcf28jiRNVOgiErTCshquenw1+8treHyu9jEPNyp0EQlKQWk1Mxd9SHlNPc9cN5XMwcleR5IjqNBF5Jj2l9cw+7HVlNfU89z10xg3sKfXkaQFMcGsZGbTzWybmeWZ2d0t3N/TzP5sZuvNbLOZzQt9VBHxwoGKWmYt+pDSKh9PXzdVZR7GjlnoZhYLPAScD4wGZpnZ6CNWuxnY4pybAJwJ/NrM4kOcVUTaWVFlLbMe+5DiyjqevDaLiYN6eR1JvkQwI/QsIM85l++c8wFLgRlHrOOARGs8oqA7cAhoCGlSEWlXJYfrmP3Yaj4tr+XJa7OYnJ7kdSQ5hmAKfQCwt9lyQdNtzT0IjAL2ARuBW51zgSP/kJnNN7NsM8suLi4+zsgi0tYOVfm46vHVFJRWs/iaKUzRBtCIEEyht3Qcrzti+TwgF+gPTAQeNLMeX/gl5xY55zKdc5mpqdrdSSQcVdTWc9Xjq9lVUsUf5k5h2tCTvI4kQQqm0AuAQc2WB9I4Em9uHrDcNcoDdgGnhCaiiLSX2no/1z+VzY6iSh69ejKnD0vxOpK0QjCFvgYYbmZDmjZ0zgRePWKdPcDZAGbWBxgJ5IcyqIi0rQZ/gP96fh0f7TrEr745gTNH9vY6krTSMfdDd841mNktwAogFljsnNtsZjc23b8Q+CnwpJltpHGK5i7nXEkb5haREHLO8T8vb+LNLQf44UWjmTHxyM1kEgmCOrDIOfc68PoRty1s9vM+4OuhjSYi7eWXK7axLHsvC742jHmnD/E6jhynoA4sEpHo9Yd/7OLhd3YyKyuNO84d4XUcOQEqdJEO7E/rCvnpX7YwfUxffnbxWF2cIsKp0EU6qH/sKOHOF9czbWgyD8ycSGyMyjzSqdBFOqDtByq56ZkcTk7tzqI5mboGaJRQoYt0MMWVdcx7Yg2d42NZPG8KPTrHeR1JQkSFLtKB1Pj8XL8km4NVdfxhbiYDenXxOpKEkM6HLtJBBAKO25flsqGgjEevmsz4gTpzYrTRCF2kg7jvja28sflTvn/haL4+pq/XcaQNqNBFOoBnV+/m0VX5zDk1nWtPH+x1HGkjKnSRKPfejmLufWUzZ41M5d7/HK19zaOYCl0kiu0qqeLmZ9cyvHd3fj87g06xestHM/3rikSpitp6bliSTWyM8dicTLonaB+IaKd/YZEo5A84bluayyclVTx93VQGJXf1OpK0A43QRaLQr97cxt+3FvHDb4zh1JN1xaGOQoUuEmVeyS3kkXd2MntqGldPS/c6jrQjFbpIFNlQUMb3/riBrCHJ/OiiMV7HkXamQheJEkUVtcxfkkNK9wQeuTKD+E56e3c02igqEgV8DQG+8+xaymvqeemm0zipe4LXkcQDKnSRKPC/r20he3cpD86exOj+PbyOIx7RZzKRCPdSTgFPfbCbG74yhP8c39/rOOIhFbpIBNtUWM49L29k2tBk7pp+itdxxGMqdJEIVVrl48ZnckjuFs+DOqxf0By6SETyBxy3LsulqKKOF248lRRtBBVU6CIR6TdvbWfV9mJ+fuk4Jg7ShSqkkT6jiUSYFZs/5cG385g5ZRCzstK8jiNhRIUuEkF2lVRx5wvrGT+wJz/6ho4ElX+nQheJEDU+Pzc9k0NsrPHwlRl0jov1OpKEGc2hi0SIe1/ZxLYDlTxxzRQGJul0uPJFGqGLRIBla/bwYk4BC84axpkje3sdR8KUCl0kzG0qLOcHr2zmK8NTuPWcEV7HkTAWVKGb2XQz22ZmeWZ291HWOdPMcs1ss5m9G9qYIh1TeU0933l2Lcld43ngionExugCz3J0x5xDN7NY4CHgXKAAWGNmrzrntjRbpxfwMDDdObfHzPSZUOQEOee488X17CurYdm3T9UZFOWYghmhZwF5zrl855wPWArMOGKd2cBy59weAOdcUWhjinQ8j67K560tB7jnglFMTk/yOo5EgGAKfQCwt9lyQdNtzY0AkszsHTPLMbM5Lf0hM5tvZtlmll1cXHx8iUU6gNX5B/nlim1cOK4f804f7HUciRDBFHpLk3buiOVOwGTgQuA84Adm9oWtN865Rc65TOdcZmpqaqvDinQEJYfrWPD8OtKSu/KLy8ZhpnlzCU4w+6EXAIOaLQ8E9rWwTolzrgqoMrNVwARge0hSinQQ/oDjtqW5lNfU8+S8LBI7x3kdSSJIMCP0NcBwMxtiZvHATODVI9Z5BfiKmXUys67AVODj0EYViX4P/j2Pf+SV8JMZY3TlIWm1Y47QnXMNZnYLsAKIBRY75zab2Y1N9y90zn1sZm8AG4AA8LhzblNbBheJNu/nlfDAyu1cOmkA38ocdOxfEDmCOXfkdHj7yMzMdNnZ2Z48tki4Kaqo5YLfvUdS13heueV0usbrrBzSMjPLcc5ltnSfXjUiHmvwB1jw/Dqq6vw8f0OGylyOm145Ih77zd+2s3rXIX5zxQSG90n0Oo5EMJ3LRcRDb28r4qG3dzIraxCXTBrodRyJcCp0EY/sK6vhjmW5nNI3kR9epItVyIlToYt4oL5p3tzXENDFKiRkNIcu4oFfrdhGzu5Sfj9rEkNTu3sdR6KERugi7Wzlxwd4dFU+V01L46IJ/b2OI1FEhS7SjgpKq7njhfWM6d+D71842us4EmVU6CLtxNcQ4Jbn1uEPOB6arXlzCT3NoYu0k/ve2Eru3jIevjKDwSndvI4jUUgjdJF2sGLzp/zhH7uYe2o6F4zr53UciVIqdJE2tudgNXe+uJ7xA3tyz4WjvI4jUUyFLtKGauv9fOe5HAx4aHYGCZ00by5tR3PoIm3oZ69tYVNhBY/NyWRQclev40iU0whdpI28klvIMx/uYf4ZQzl3dB+v40gHoEIXaQN5RYf57+UbyUxP4rvnjfQ6jnQQKnSREKvx+bn52bV0jovl97MnERert5m0D82hi4TYD17ZxPaiSp6al0W/nl28jiMdiIYOIiH0QvZe/phTwIKzhnHGiFSv40gHo0IXCZHN+8r5wZ82cdrJJ3HrOSO8jiMdkApdJATKa+q56Zm19Ooax+9mTSI2xryOJB2Q5tBFTpBzjjtfXM++shqWfXsaKd0TvI4kHZRG6CIn6NFV+by15QD3XDCKyenJXseRDkyFLnICPth5kPvf2MqF4/sx7/TBXseRDk6FLnKciipqWfD8OgandOO+y8Zjpnlz8Zbm0EWOQ70/wM3PraWqroHnbphK9wS9lcR7ehWKHIf739jKmk9K+e3MiYzok+h1HBFAUy4irfbn9ft47L1dzDk1nRkTB3gdR+RzKnSRVtj6aQXf++MGMtOTdJFnCTsqdJEglVfX8+2nc0js3ImHr8wgvpPePhJeNIcuEoRAwHHbsnXsK6th6fxp9O7R2etIIl8Q1BDDzKab2TYzyzOzu79kvSlm5jezy0MXUcR7D6zcwdvbirn3ojE6eEjC1jEL3cxigYeA84HRwCwz+8LkYdN69wErQh1SxEtvbTnA71bu4PLJA7lqaprXcUSOKpgRehaQ55zLd875gKXAjBbWWwC8BBSFMJ+Ip3YWH+aOZbmMG9CTn108VgcPSVgLptAHAHubLRc03fY5MxsAXAIs/LI/ZGbzzSzbzLKLi4tbm1WkXVXWNm4EjesUw8KrJ9M5LtbrSCJfKphCb2lI4o5YfgC4yznn/7I/5Jxb5JzLdM5lpqbq5P8SvvwBx21Lc9lVUsWDsyYxoJeuPCThL5i9XAqAQc2WBwL7jlgnE1ja9HE0BbjAzBqcc38KSUqRdvbrN7excmsRP5kxhtOGpXgdRyQowRT6GmC4mQ0BCoGZwOzmKzjnhnz2s5k9CfxFZS6R6pXcQh5+ZyezstK4elq613FEgnbMQnfONZjZLTTuvRILLHbObTazG5vu/9J5c5FIsn5vGd/74wayBifz42+M0UZQiShBHVjknHsdeP2I21oscufcNSceS6T9FVXUMv/pbFK6J/DIVToSVCKPjhQVAWrr/cx/OofK2gZeuuk0TtJl5CQCqdClw3POcc/yjeTuLWPhVRmM6tfD60gix0WfKaXDe+TdnSxfV8jt54xg+th+XscROW4qdOnQ/rJhH/e/sY1vTOjPgq8N8zqOyAlRoUuHlbO7lDteWE9mehL3Xz6emBjt0SKRTYUuHdLug1XcsCSb/j07s2hOpg7rl6igQpcOp6zax7wn1xBwjifmZZHcLd7rSCIhoUKXDsXXEODbT+dQcKiGRVdnMiSlm9eRREJGuy1Kh+Gc4+7lG1i96xAPXDGRrCG6UIVEF43QpcP4zd92sHxt4+6JF08acOxfEIkwKnTpEJ7+cDe/W7mDb04eyH+drd0TJTqp0CXqvb5xP/e+somzT+nNzy8dpxNuSdRSoUtU+2DnQW5bmktGWhIPzs6gU6xe8hK99OqWqLV5Xznzl2QzOKUrf5ibSZd47Wsu0U2FLlFpz8Fq5i5eQ2LnTjx1bRa9umpfc4l+KnSJOsWVdVy9eDUNgQBLrsuiX09dD1Q6BhW6RJXy6nrmLv6IAxW1LL5mCsN6J3odSaTdqNAlalTW1jPniY/IKzrMo1dnkpGW5HUkkXalQpeoUFXXwLwn1rC5sJyHr8zgqyNSvY4k0u5U6BLxauv9XP9UNmv3lPLbmZM4Z3QfryOJeELncpGIVtfg59tP5/DhroP837cmcOF4XXFIOi6N0CVi1fsD3PLcOt7dXszPLxnHJZMGeh1JxFMqdIlI9f4Aty3N5a0tB/jxN8YwMyvN60gintOUi0ScugY/C55bx5tbDvA/F4xi7mmDvY4kEhZU6BJRauv93PhMDu9sK+ZHF43mmtOHeB1JJGyo0CViVPsauP6pbD7IP8jPLx3HLE2ziPwbFbpEhMraeuY9sYa1e0r59TcncGmGNoCKHEmFLmGvvLrxCNDNheX8flaGdk0UOQoVuoS14so65i5uPJz/kasmc64OGhI5KhW6hK384sPMfeIjSip9PDY3U4fzixxDUPuhm9l0M9tmZnlmdncL919pZhuavv5pZhNCH1U6knV7Srl84QdU1fl5fv40lblIEI45QjezWOAh4FygAFhjZq8657Y0W20X8FXnXKmZnQ8sAqa2RWCJfis/PsDNz62ld2Jnnro2iyEp3byOJBIRghmhZwF5zrl855wPWArMaL6Cc+6fzrnSpsUPAe2CIMdl6Ud7uGFJNsN7J/LSTaepzEVaIZg59AHA3mbLBXz56Ps64K8t3WFm84H5AGlp2odY/sU5x29X7uCBv+3gjBGpPHJlBt0StIlHpDWCecdYC7e5Flc0O4vGQv+Plu53zi2icTqGzMzMFv+GdDy19X7+e/lGXl5XyGUZA/nFZeOIi9VphkRaK5hCLwAGNVseCOw7ciUzGw88DpzvnDsYmngS7Q5U1DL/6RzW7y3jjnNHsOBrwzBraQwhIscSTKGvAYab2RCgEJgJzG6+gpmlAcuBq51z20OeUqJS7t4y5i/J5nBdAwuvmsz0sX29jiQS0Y5Z6M65BjO7BVgBxAKLnXObzezGpvsXAvcCJwEPN42uGpxzmW0XWyLd8rUF3L18I70TE1h+3Wmc0reH15FEIp45581UdmZmpsvOzvbkscU7/oDjvje2smhVPtOGJvPwlZNJ7hbvdSyRiGFmOUcbMGs3Amk3RZW13L4sl/fzDnL1tHTuvWi0Nn6KhJAKXdrF+3kl3Lo0l8raeu67bBxXTNFuqyKhpkKXNuUPNO5f/vu/72BoSjeevX4qI/smeh1LJCqp0KXNHKio5dal6/gw/xCXZQzkpxePoWu8XnIibUXvLmkT72wr4v+9sJ5qn59ffXMCl0/W2SBE2poKXULqcF0D//vaxzz/0R5G9OnO0tkZDO+jKRaR9qBCl5D5584SvvfHDRSW1fDtM4Zy+7kj6BwX63UskQ5DhS4nrNrXwP1vbOPJf37CkJRu/PHGU5mcnux1LJEOR4UuJ2TNJ4f47ovr+eRgNdecNpi7pp9Cl3iNykW8oEKX43Koysd9f93Ksuy9DEruwtL505g29CSvY4l0aCp0aZVAwLF0zV7uX7GVw7UNzD9jKLeePVznLhcJA3oXStA2FpTz/Vc2sX5vGVlDkvnZxWMZoT1YRMKGCl2O6VCVj9+8tZ1nVu/mpG4J/OaKCVw8cYDOWy4SZlToclTVvgYW/2MXC9/Np9rXwJxp6dzx9ZH07BLndTQRaYEKXb6gwR9gWfZeHvjbDoor6zh3dB/umj6SYb01vSISzlTo8jnnHCs2f8r9K7aRX1zF5PQkHrkyg8zB2qdcJBKo0AV/wPHXTft56O2dfLy/gmG9u/PYnEzOGdVb8+QiEUSF3oHV+wO8vK6Qhe/sJL+kiqGp3fjVNydw8cT+dNKFJ0Qijgq9A6rx+Xkhey+LVuVTWFbD6H49ePjKDM4b05fYGI3IRSKVCr0D2XOwmmdW72bZmr2U19STmZ7Ezy4Zy5kjUjW1IhIFVOhRLhBwvLujmKc/2M3b24qIMeO8MX2Ye+pgsoYkq8hFoogKPUodqKjlT+sKef6jPXxysJqU7gks+NpwZmel0bdnZ6/jiUgbUKFHkWpfA29uPsBLawt4P6+EgIPM9CTu+PpIpo/pS3wnbegUiWYq9AhX7w+wOv8QL68r5I1N+6ny+RnQqws3nzWMSyYNYGhqd68jikg7UaFHoBqfn1U7ilmx+VNWflxEeU09iQmd+M/x/bk0YwBTBicTo71VRDocFXqEKK6sY9X2Yt7c8invbi+mtj5Aj86dOGdUH74+pg9njuyty72JdHAq9DBV4/OzetdB3s8r4b0dJWz9tBKAvj06863MQZw3pi9ZQ5KJ0wFAItJEhR4mKmvryd1bRs7uUlbnHyJndyk+f4D42BgyByfxvekj+cqwVMb076HpFBFpkQrdA4GAY9fBKjYUlJH9SSk5u0vZdqAS58AMTunbg7mnpfMfw1PJGpysa3SKSFBU6G2stt7P9gOVbN5XwZZ9FWzeV87WTyup9vkB6J7QiUlpvZg+ti+T05OYOKgXiZ11vnERaT0VeggEAo6iyjrySw6TX1zFzuJ/fS8sq8G5xvUSEzoxqn8PvpU5iNH9ezBuQE9G9EnU+VNEJCSCKnQzmw78FogFHnfO/eKI+63p/guAauAa59zaEGf1RCDgKKupp7iyjuLKOg5U1FJYVkNhaQ0FZdUUltawr6wWnz/w+e90jY9laGo3MtKSuHzyQEb2SWRM/54MTOqi+W8RaTPHLHQziwUeAs4FCoA1Zvaqc25Ls9XOB4Y3fU0FHmn67innHP6Aw+cPUOPzU+3zU+VroKrOT7WvgWqfn8raBsqqfZTX1FNWXd/4vaae0iofxZV1lByuoyHgvvC3UxMTGJjUhbEDejJ9bD8GJHVhyEndOLl3N/r26KxzpIhIuwtmhJ4F5Dnn8gHMbCkwA2he6DOAJc45B3xoZr3MrJ9zbn+oA7+zrYifvfYxgYAj4BwBBwHncK7xQg31/gA+fwBfQ+N398UublGMQY8ucfTqEkfPrvEkd4vnlL6JpCYmkJqYQO/Ezp//3K9nZ+3zLSJhJ5hCHwDsbbZcwBdH3y2tMwD4t0I3s/nAfIC0tLTWZgUgsXMcI/skYgaxMUaMGWYQY0aMQVxsDHGxMSR0iiG+U8zny13jY+kaH0u3hE7/9j0xIY6eXeNITOik6RARiWjBFHpLLXfkuDeYdXDOLQIWAWRmZgY5dv53k9OTmJyedDy/KiIS1YI5zLAAGNRseSCw7zjWERGRNhRMoa8BhpvZEDOLB2YCrx6xzqvAHGs0DShvi/lzERE5umNOuTjnGszsFmAFjbstLnbObTazG5vuXwi8TuMui3k07rY4r+0ii4hIS4LaD9059zqNpd38toXNfnbAzaGNJiIiraFT9YmIRAkVuohIlFChi4hECRW6iEiUMBfssfGhfmCzYmD3cf56ClASwjihEq65IHyzKVfrKFfrRGOudOdcakt3eFboJ8LMsp1zmV7nOFK45oLwzaZcraNcrdPRcmnKRUQkSqjQRUSiRKQW+iKvAxxFuOaC8M2mXK2jXK3ToXJF5By6iIh8UaSO0EVE5AgqdBGRKBHWhW5m081sm5nlmdndLdx/ipl9YGZ1ZnZnGOW60sw2NH3908wmhEmuGU2Zcs0s28z+IxxyNVtvipn5zezycMhlZmeaWXnT85VrZveGQ65m2XLNbLOZvRsOuczsu82eq01N/5bJYZCrp5n92czWNz1f7XI22CByJZnZy03vyY/MbOwJP6hzLiy/aDxV705gKBAPrAdGH7FOb2AK8L/AnWGU6zQgqenn84HVYZKrO//abjIe2BoOuZqt93caz+p5eTjkAs4E/tIer6tW5upF4zV905qWe4dDriPWvwj4ezjkAu4B7mv6ORU4BMSHQa5fAj9s+vkUYOWJPm44j9A/vzi1c84HfHZx6s8554qcc2uA+jDL9U/nXGnT4oc0XsEpHHIddk2vHqAbLVwm0ItcTRYALwFF7ZCpNbnaWzC5ZgPLnXN7oPF9ECa5mpsFPB8muRyQaGZG46DmENAQBrlGAysBnHNbgcFm1udEHjScC/1oF572WmtzXQf8tU0TNQoql5ldYmZbgdeAa8Mhl5kNAC4BFtJ+gv13PLXpo/pfzWxMmOQaASSZ2TtmlmNmc8IkFwBm1hWYTuP/oMMh14PAKBovi7kRuNU5FwiDXOuBSwHMLAtI5wQHf+Fc6EFdeNoDQecys7NoLPS72jRR08O1cFtLF+p+2Tl3CnAx8NM2TxVcrgeAu5xz/nbI85lgcq2l8bwZEwXYL44AAAH9SURBVIDfA39q81TB5eoETAYuBM4DfmBmI8Ig12cuAt53zh1qwzyfCSbXeUAu0B+YCDxoZj3CINcvaPwfcy6Nn1DXcYKfHIK6YpFHwvXC00HlMrPxwOPA+c65g+GS6zPOuVVmdrKZpTjn2vLkRcHkygSWNn4iJgW4wMwanHNtWaDHzOWcq2j28+tm9nCYPF8FQIlzrgqoMrNVwARgu8e5PjOT9plugeByzQN+0TTdmGdmu2ics/7Iy1xNr695AE3TQbuavo5fW2+0OIGNCp2AfGAI/9qoMOYo6/6I9tsoesxcQBqN11c9LZyeL2AY/9oomgEUfrYcDv+OTes/SftsFA3m+erb7PnKAvaEw/NF4/TByqZ1uwKbgLFe52paryeNc9Td2vrfsBXP1yPAj5p+7tP0uk8Jg1y9aNo4C9wALDnRxw3bEboL4uLUZtYXyAZ6AAEzu43GLckVR/3D7ZALuBc4CXi4adTZ4Nr4jG9B5roMmGNm9UANcIVrejV5nKvdBZnrcuAmM2ug8fmaGQ7Pl3PuYzN7A9gABIDHnXObvM7VtOolwJuu8dNDmwsy10+BJ81sI41TIXe5tv2UFWyuUcASM/PTuNfSdSf6uDr0X0QkSoTzRlEREWkFFbqISJRQoYuIRAkVuohIlFChi4hECRW6iEiUUKGLiESJ/w8hvp1TJ/0RGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "krw_plot = [krw(sw) for sw in sw_plot]\n",
    "kro_plot = [kro(sw) for sw in sw_plot]\n",
    "fw_plot = [fw(sw) for sw in sw_plot]\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(sw_plot, krw_plot, sw_plot, kro_plot)\n",
    "plt.show()\n",
    "\n",
    "plt.figure(2)\n",
    "plt.plot(sw_plot, fw_plot)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create the grid\n",
    "mesh = Grid1D(dx = Lx/nx, nx = nx)\n",
    "x = mesh.cellCenters\n",
    "\n",
    "# create the cell variables and boundary conditions\n",
    "sw = CellVariable(mesh=mesh, name=\"saturation\", hasOld=True, value = swc)\n",
    "p = CellVariable(mesh=mesh, name=\"pressure\", hasOld=True, value = p0)\n",
    "# sw.setValue(1,where = x<=dx)\n",
    "sw.constrain(1.0,mesh.facesLeft)\n",
    "#sw.constrain(0., mesh.facesRight)\n",
    "sw.faceGrad.constrain([0], mesh.facesRight)\n",
    "p.constrain(p0, mesh.facesRight)\n",
    "p.constrain(1.02*p0, mesh.facesLeft)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equations\n",
    "$$\\varphi \\frac{\\partial S_w}{\\partial t}+u \\frac{\\partial f_w}{\\partial x}=0$$ or\n",
    "$$\\varphi \\frac{\\partial S_w}{\\partial t}+\\nabla.\\left( u \\frac{\\partial f_w}{\\partial S_w} S_w\\right)+ \\nabla. \\left( u f_w-u\\frac{\\partial f_w}{\\partial S_w} S_{w0} \\right)=0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0005369273839403227\n",
      "0.0005369273839403227\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/ipykernel_launcher.py:22: RuntimeWarning: invalid value encountered in greater\n",
      "/home/ali/miniconda3/envs/myfipy/lib/python3.7/site-packages/ipykernel_launcher.py:23: RuntimeWarning: invalid value encountered in less\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "Factor is exactly singular",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-113-986a9bfede1b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     19\u001b[0m     \u001b[0mloop_count\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m     \u001b[0;32mwhile\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m         \u001b[0mswres_new\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0meq\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msweep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     22\u001b[0m         \u001b[0msw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m>\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0msor\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0msor\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m         \u001b[0msw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m<\u001b[0m\u001b[0mswc\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mswc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/terms/term.py\u001b[0m in \u001b[0;36msweep\u001b[0;34m(self, var, solver, boundaryConditions, dt, underRelaxation, residualFn, cacheResidual, cacheError)\u001b[0m\n\u001b[1;32m    230\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mresidualVector\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    231\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 232\u001b[0;31m         \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_solve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    233\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    234\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mresidual\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/solvers/scipy/scipySolver.py\u001b[0m in \u001b[0;36m_solve\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m     24\u001b[0m              \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"SciPy solvers cannot be used with multiple processors\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     25\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m          \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnumerix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_solve_\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mravel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnumerix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mRHSvector\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/miniconda3/envs/myfipy/lib/python3.7/site-packages/fipy/solvers/scipy/linearLUSolver.py\u001b[0m in \u001b[0;36m_solve_\u001b[0;34m(self, L, x, b)\u001b[0m\n\u001b[1;32m     32\u001b[0m                                             \u001b[0mrelax\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     33\u001b[0m                                             \u001b[0mpanel_size\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 34\u001b[0;31m                                             permc_spec=3)\n\u001b[0m\u001b[1;32m     35\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     36\u001b[0m         \u001b[0merror0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnumerix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnumerix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mL\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mx\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/myfipy/lib/python3.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py\u001b[0m in \u001b[0;36msplu\u001b[0;34m(A, permc_spec, diag_pivot_thresh, relax, panel_size, options)\u001b[0m\n\u001b[1;32m    309\u001b[0m         \u001b[0m_options\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    310\u001b[0m     return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,\n\u001b[0;32m--> 311\u001b[0;31m                           ilu=False, options=_options)\n\u001b[0m\u001b[1;32m    312\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    313\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mRuntimeError\u001b[0m: Factor is exactly singular"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAEICAYAAABWCOFPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAT5UlEQVR4nO3dfZBU5Z3F8e8RkBdFjczogmAGa9GIvGcQzEYXdUUQI27FbMAXgrFC4cuapGIibmqFmFrLVLkxsSSybNZgXBU3aBGCrEYjbCRq1p5EUUCUKOqI0QE3iCilyG//6Dukn3aGaaCnu4nnUzXF3HufvvdMQ595+nbTVxGBmVmrA6odwMxqi0vBzBIuBTNLuBTMLOFSMLOES8HMEi4FqxhJF0j6ZbVz2O7J71P4eJIUwKCIWN9J+28AXgK6RcSOzjiGdQ7PFGyvSOpS7QzWOVwKfwEkXS3pNUlbJa2TdLqkEyU9LulPkl6XdIukA7Pxv85u+rSkdyR9UdJ0SSuL9huS/jr7foGkWyUtk7QNOFXSJEm/l/S2pFclzSm4eesx/pQd46TiY0j6jKQnJW3J/vxMwbYVkr4r6TfZz/VLSXWdcPdZEZfCfk7SccAVwOiI6A2cCWwAPgS+DtQBJwGnA5cBRMQp2c2HR8TBEXFPiYc7H/gXoDewEtgGTAMOAyYBl0o6NxvbeozDsmM8XpT7cOB+4GagD/B94H5JfYqOdzFwBHAgcFWJOW0fuBT2fx8C3YHBkrpFxIaI+ENENEXEExGxIyI2AP8G/O0+HuvnEfGbiNgZEdsjYkVEPJMtrwLu3oNjTAJeiIg7sox3A88BnysY85OIeD4i3gP+Cxixj/mtBC6F/Vx2ovBrwBzgTUkLJfWTdKykpZL+KOlt4Hrys4Z98WrhgqQxkpZLapG0BZi5B8foB7xctO5l4KiC5T8WfP8ucPAe5rW94FL4CxARd0XEZ4FPAgF8D7iV/G/eQRFxCPBPgHazm21Ar9YFSX/V1qGKlu8ClgADIuJQYF7BMTp6WWtjlrfQ0cBrHdzOOplLYT8n6ThJp0nqDmwH3iP/lKI38DbwjqRPAZcW3fQN4JiC5aeBEySNkNSD/MyjI72BtyJiu6QTyZ8DaNUC7Cw6RqFlwLGSzpfUVdIXgcHA0hKOa53IpbD/6w7cAGwiP90+gvys4CryD9KtwL8DxScT5wC3Z69O/ENEPA9cBzwMvED+RGJHLgOuk7QVuJb8834AIuJd8iclf5MdY2zhDSNiM3A28A1gM/At4OyI2FT6j26dwW9eMrOEZwpmlnApmFnCpWBmCZeCmSW6VuvAdXV10dDQUK3Dm31sNTU1bYqI+va2V60UGhoayOVy1Tq82ceWpOJ3kib89MHMEi4FM0u4FMwsUbVzCm354IMPaG5uZvv27dWO8helR48e9O/fn27dulU7iu0HaqoUmpub6d27Nw0NDUi7+w99VqqIYPPmzTQ3NzNw4MBqx7H9QE09fdi+fTt9+vRxIZSRJPr06ePZl5WspkoBcCF0At+ntidqrhTMrLpcCvtgwYIFbNy4sWz727BhA3fdddeu5Vwux5VXXlm2/ZuVwqWwD/amFHbsaP+6KMWl0NjYyM0337zX+cz2hkuhyLZt25g0aRLDhw9nyJAh3HPPPVx33XWMHj2aIUOGMGPGDCKCRYsWkcvluOCCCxgxYgTvvfceDQ0NbNqU/+CgXC7HuHHjAJgzZw4zZsxg/PjxTJs2jQ0bNnDyySczatQoRo0axWOPPQbArFmzePTRRxkxYgQ33XQTK1as4Oyzzwbgrbfe4txzz2XYsGGMHTuWVatW7dr3l7/8ZcaNG8cxxxzjErF9VlMvSRb6zi9Ws2bj22Xd5+B+hzD7cyfsdswDDzxAv379uP/++wHYsmULZ5xxBtdeey0AF110EUuXLuW8887jlltu4cYbb6SxsbHDYzc1NbFy5Up69uzJu+++y0MPPUSPHj144YUXmDp1KrlcjhtuuIEbb7yRpUvzH1O4YsWKXbefPXs2I0eOZPHixTzyyCNMmzaNp556CoDnnnuO5cuXs3XrVo477jguvfRSvyfB9ppnCkWGDh3Kww8/zNVXX82jjz7KoYceyvLlyxkzZgxDhw7lkUceYfXq1Xu833POOYeePXsC+TdpfeUrX2Ho0KF84QtfYM2aNR3efuXKlVx00UUAnHbaaWzevJktW7YAMGnSJLp3705dXR1HHHEEb7zxxh7nM2vV4UxB0m3kP2DzzYgY0sZ2AT8EziL/2fzTI+J3+xqso9/oneXYY4+lqamJZcuWcc011zB+/Hjmzp1LLpdjwIABzJkzp93X/Lt27crOnTsBPjLmoIMO2vX9TTfdxJFHHsnTTz/Nzp076dGjR4e52voszdaXGrt3775rXZcuXXZ73sKsI6XMFBYAE3azfSIwKPuaQf56A/utjRs30qtXLy688EKuuuoqfve7fL/V1dXxzjvvsGjRol1je/fuzdatW3ctNzQ00NTUBMC9997b7jG2bNlC3759OeCAA7jjjjv48MMP29xfoVNOOYU777wTyD+tqKur45BDDtm3H9asDR3OFCLi19llxdszGfhp5H+VPSHpMEl9I+L1MmWsqGeeeYZvfvObHHDAAXTr1o1bb72VxYsXM3ToUBoaGhg9evSusdOnT2fmzJn07NmTxx9/nNmzZ3PJJZdw/fXXM2bMmHaPcdlll/H5z3+en/3sZ5x66qm7ZhHDhg2ja9euDB8+nOnTpzNy5Mhdt5kzZw4XX3wxw4YNo1evXtx+++2ddyfYx1pJH/GelcLSdp4+LAVuiIiV2fKvgKsj4iOfoCJpBvnZBEcfffSnX345/ayHtWvXcvzxx+/5T2Ed8n1rrSQ1RUS7Z8fLcaKxrffQttk0ETE/IhojorG+vt1PgzKzKipHKTQDAwqW+5O/TqCZ7YfKUQpLgGnKGwts2ZfzCb5iVfn5PrU9UcpLkncD44A6Sc3AbKAbQETMI3+h0LOA9eRfkrx4b8P06NGDzZs3+79Pl1Hr5ymU8rKnGZT26sPUDrYHcHk5wvTv35/m5mZaWlrKsTvLtH7yklkpauptzt26dfOnA5lVmd/mbGYJl4KZJVwKZpZwKZhZwqVgZgmXgpklXApmlnApmFnCpWBmCZeCmSVcCmaWcCmYWcKlYGYJl4KZJVwKZpZwKZhZwqVgZgmXgpklXApmlnApmFnCpWBmCZeCmSVcCmaWcCmYWcKlYGYJl4KZJVwKZpYoqRQkTZC0TtJ6SbPa2H6opF9IelrSakl7feVpM6uuDktBUhdgLjARGAxMlTS4aNjlwJqIGE7+svX/KunAMmc1swooZaZwIrA+Il6MiPeBhcDkojEB9JYk4GDgLWBHWZOaWUWUUgpHAa8WLDdn6wrdAhwPbASeAb4aETuLdyRphqScpFxLS8teRjazzlRKKaiNdVG0fCbwFNAPGAHcIumQj9woYn5ENEZEY319/R6HNbPOV0opNAMDCpb7k58RFLoYuC/y1gMvAZ8qT0Qzq6RSSuFJYJCkgdnJwynAkqIxrwCnA0g6EjgOeLGcQc2sMrp2NCAidki6AngQ6ALcFhGrJc3Mts8DvgsskPQM+acbV0fEpk7MbWadpMNSAIiIZcCyonXzCr7fCIwvbzQzqwa/o9HMEi4FM0u4FMws4VIws4RLwcwSLgUzS7gUzCzhUjCzhEvBzBIuBTNLuBTMLOFSMLOES8HMEi4FM0u4FMws4VIws4RLwcwSLgUzS7gUzCzhUjCzhEvBzBIuBTNLuBTMLOFSMLOES8HMEi4FM0u4FMwsUVIpSJogaZ2k9ZJmtTNmnKSnJK2W9D/ljWlmldLhBWYldQHmAmcAzcCTkpZExJqCMYcBPwImRMQrko7orMBm1rlKmSmcCKyPiBcj4n1gITC5aMz5wH0R8QpARLxZ3phmVimllMJRwKsFy83ZukLHAp+QtEJSk6Rpbe1I0gxJOUm5lpaWvUtsZp2qlFJQG+uiaLkr8GlgEnAm8M+Sjv3IjSLmR0RjRDTW19fvcVgz63wdnlMgPzMYULDcH9jYxphNEbEN2Cbp18Bw4PmypDSziillpvAkMEjSQEkHAlOAJUVjfg6cLKmrpF7AGGBteaOaWSV0OFOIiB2SrgAeBLoAt0XEakkzs+3zImKtpAeAVcBO4McR8WxnBjezzqGI4tMDldHY2Bi5XK4qxzb7OJPUFBGN7W33OxrNLOFSMLOES8HMEi4FM0u4FMws4VIws4RLwcwSLgUzS7gUzCzhUjCzhEvBzBIuBTNLuBTMLOFSMLOES8HMEi4FM0u4FMws4VIws4RLwcwSLgUzS7gUzCzhUjCzhEvBzBIuBTNLuBTMLOFSMLOES8HMEiWVgqQJktZJWi9p1m7GjZb0oaTzyhfRzCqpw1KQ1AWYC0wEBgNTJQ1uZ9z3yF+d2sz2U6XMFE4E1kfEixHxPrAQmNzGuH8E7gXeLGM+M6uwUkrhKODVguXmbN0uko4C/h6Yt7sdSZohKScp19LSsqdZzawCSikFtbEuipZ/AFwdER/ubkcRMT8iGiOisb6+vtSMZlZBXUsY0wwMKFjuD2wsGtMILJQEUAecJWlHRCwuS0ozq5hSSuFJYJCkgcBrwBTg/MIBETGw9XtJC4ClLgSz/VOHpRAROyRdQf5VhS7AbRGxWtLMbPtuzyOY2f6llJkCEbEMWFa0rs0yiIjp+x7LzKrF72g0s4RLwcwSLgUzS7gUzCzhUjCzhEvBzBIuBTNLuBTMLOFSMLOES8HMEi4FM0u4FMws4VIws4RLwcwSLgUzS7gUzCzhUjCzhEvBzBIuBTNLuBTMLOFSMLOES8HMEi4FM0u4FMws4VIws4RLwcwSLgUzS5RUCpImSFonab2kWW1sv0DSquzrMUnDyx/VzCqhw1KQ1AWYC0wEBgNTJQ0uGvYS8LcRMQz4LjC/3EHNrDJKmSmcCKyPiBcj4n1gITC5cEBEPBYR/5ctPgH0L29MM6uUUkrhKODVguXmbF17LgH+u60NkmZIyknKtbS0lJ7SzCqmlFJQG+uizYHSqeRL4eq2tkfE/IhojIjG+vr60lOaWcV0LWFMMzCgYLk/sLF4kKRhwI+BiRGxuTzxzKzSSpkpPAkMkjRQ0oHAFGBJ4QBJRwP3ARdFxPPlj2lmldLhTCEidki6AngQ6ALcFhGrJc3Mts8DrgX6AD+SBLAjIho7L7aZdRZFtHl6oNM1NjZGLperyrHNPs4kNe3ul7bf0WhmCZeCmSVcCmaWcCmYWcKlYGYJl4KZJVwKZpZwKZhZwqVgZgmXgpklXApmlnApmFmilM9TqJrv/GI1aza+Xe0YZvuFwf0OYfbnTtjn/XimYGaJmp4plKP1zGzPeKZgZgmXgpklXApmlnApmFnCpWBmCZeCmSVcCmaWcCmYWcKlYGYJl4KZJVwKZpZwKZhZwqVgZomSSkHSBEnrJK2XNKuN7ZJ0c7Z9laRR5Y9qZpXQYSlI6gLMBSYCg4GpkgYXDZsIDMq+ZgC3ljmnmVVIKTOFE4H1EfFiRLwPLAQmF42ZDPw08p4ADpPUt8xZzawCSimFo4BXC5abs3V7OgZJMyTlJOVaWlr2NKuZVUAppaA21sVejCEi5kdEY0Q01tfXl5LPzCqslFJoBgYULPcHNu7FGDPbD5RSCk8CgyQNlHQgMAVYUjRmCTAtexViLLAlIl4vc1Yzq4AOP7g1InZIugJ4EOgC3BYRqyXNzLbPA5YBZwHrgXeBizsvspl1ppI+zTkilpF/4Beum1fwfQCXlzeamVWD39FoZgmXgpklXApmlnApmFnCpWBmCZeCmSVcCmaWcCmYWcKlYGYJl4KZJVwKZpZwKZhZQvn/y1SFA0stwMsdDKsDNlUgzr5wxvJwxvIoJeMnI6LdTzmqWimUQlIuIhqrnWN3nLE8nLE8ypHRTx/MLOFSMLNErZfC/GoHKIEzloczlsc+Z6zpcwpmVnm1PlMwswpzKZhZomZLoaOL2laDpAGSlktaK2m1pK9m6w+X9JCkF7I/P1HlnF0k/V7S0lrMl2U6TNIiSc9l9+dJtZZT0tezv+dnJd0tqUe1M0q6TdKbkp4tWNduJknXZI+hdZLOLOUYNVkKJV7Uthp2AN+IiOOBscDlWa5ZwK8iYhDwq2y5mr4KrC1YrrV8AD8EHoiITwHDyeetmZySjgKuBBojYgj5yxtMqYGMC4AJRevazJT925wCnJDd5kfZY2v3IqLmvoCTgAcLlq8Brql2rjZy/hw4A1gH9M3W9QXWVTFT/+wfxmnA0mxdzeTLMhwCvER2ortgfc3k5M/XRz2c/KUQlgLjayEj0AA829H9Vvy4IX/tlpM62n9NzhQo8YK11SSpARgJ/BY4MrIrYmV/HlG9ZPwA+Baws2BdLeUDOAZoAX6SPc35saSDqKGcEfEacCPwCvA6+aue/bKWMhZoL9NePY5qtRRKumBttUg6GLgX+FpEvF3tPK0knQ28GRFN1c7Sga7AKODWiBgJbKM2ntLskj0vnwwMBPoBB0m6sLqp9thePY5qtRRq9oK1krqRL4Q7I+K+bPUbkvpm2/sCb1Yp3t8A50jaACwETpP0nzWUr1Uz0BwRv82WF5EviVrK+XfASxHREhEfAPcBn6mxjK3ay7RXj6NaLYVSLmpbcZIE/AewNiK+X7BpCfCl7PsvkT/XUHERcU1E9I+IBvL32SMRcWGt5GsVEX8EXpV0XLbqdGANtZXzFWCspF7Z3/vp5E+G1lLGVu1lWgJMkdRd0kBgEPC/He6tWidySjiZchbwPPAH4NvVzpNl+iz56dcq4Kns6yygD/mTey9kfx5eA1nH8ecTjbWYbwSQy+7LxcAnai0n8B3gOeBZ4A6ge7UzAneTP8fxAfmZwCW7ywR8O3sMrQMmlnIMv83ZzBK1+vTBzKrEpWBmCZeCmSVcCmaWcCmYWcKlYGYJl4KZJf4fpePsr7ywunMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "eq_p = DiffusionTerm(var=p, coeff=-k*(krw(sw.faceValue)/muw+kro(sw.faceValue)/muo))- \\\n",
    "UpwindConvectionTerm(var=sw, coeff=-k*(dkrw(sw.faceValue)/muw+dkro(sw.faceValue)/muo)*p.faceGrad)- \\\n",
    "(k*(dkrw(sw.faceValue)/muw+dkro(sw.faceValue)/muo)*sw.faceValue*p.faceGrad).divergence == 0\n",
    "\n",
    "eq_sw = TransientTerm(coeff=phi, var=sw) + \\\n",
    "DiffusionTerm(var=p, coeff=-k*krw(sw.faceValue)/muw)+ \\\n",
    "UpwindConvectionTerm(var=sw, coeff=-k*dkrw(sw.faceValue)/muw*p.faceGrad)- \\\n",
    "(-k*dkrw(sw.faceValue)/muw*p.faceGrad*sw.faceValue).divergence == 0\n",
    "\n",
    "eq = eq_p & eq_sw\n",
    "steps = 200\n",
    "dt0 = 500.\n",
    "dt = dt0\n",
    "t_end = steps*dt0\n",
    "t = 0.0\n",
    "viewer = Viewer(vars = sw, datamax=1.1, datamin=-0.1)\n",
    "while t<t_end:\n",
    "    swres = 1.0e6\n",
    "    loop_count = 0\n",
    "    while True:\n",
    "        swres_new = eq.sweep(dt = dt)\n",
    "        sw.value[sw.value>1-sor]=1-sor\n",
    "        sw.value[sw.value<swc]=swc\n",
    "        loop_count+=1\n",
    "        if loop_count==1:\n",
    "            sw_res = swres_new\n",
    "        if swres_new>sw_res or loop_count>5:\n",
    "            dt = dt/3\n",
    "            sw.value = sw.old[:]\n",
    "            p.value = p.old[:]\n",
    "            continue\n",
    "        swres=swres_new\n",
    "        print(swres)\n",
    "        if swres_new<1e-5:\n",
    "            sw.updateOld()\n",
    "            p.updateOld()\n",
    "            t+=dt\n",
    "            dt = dt0\n",
    "            break\n",
    "        \n",
    "# Note: try to use the Appleyard method; the overflow is a result of wrong rel-perm values        \n",
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytical solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fractional_flow as ff\n",
    "xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj, R_oil = ff.frac_flow_wf(muw=muw, muo=muo, ut=u, phi=1.0, \\\n",
    "  k=1e-12, swc=swc, sor=sor, kro0=kro0, no=no, krw0=krw0, \\\n",
    "  nw=nw, sw0=swc, sw_inj=1.0, L=Lx, pv_inj=5.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(xt_prf, sw_prf)\n",
    "plt.plot(x.value.squeeze()/(steps*dt), sw.value)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "?sw.old"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sw.value[sw.value>1-sor]=1-sor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sw.value)\n",
    "sw.value=sw.old[:]\n",
    "print(sw.value)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
